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A DISCONTINUOUS GALERKIN FINITE ELEMENT METHOD FOR 
HAMILTON-JACOBI EQUATIONS 

CHANGQING HU AND CHI- WANG SHU * 


Abstract. In this paper, we present a discontinuous Galerkin finite element method for solving the 
nonlinear Hamilton- Jacobi equations. This method is based on the Runge-Kutta discontinuous Galerkin 
finite element method for solving conservation laws. The method has the flexibility of treating complicated 
geometry by using arbitrary triangulation, can achieve high order accuracy with a local, compact stencil, 
and are suited for efficient parallel implementation. One and two dimensional numerical examples are given 
to illustrate the capability of the method. 
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1. Introduction. In this paper we consider the numerical solutions of the Hamilton- Jacobi (HJ) equa- 
tions 

(1.1) (f> t + H{4> Xl ,...,<t> Xd ) = 0, <j)(x,0) = 4>°(x). 

As is well known, the solutions to (1.1) are Lipschitz continuous but may have discontinuous derivatives, 
regardless of the smoothness of the initial condition 4>°(x). The non- uniqueness of such solutions also 
necessitates the definition of viscosity solutions, to single out a unique, practically relevant solution. (See 
Crandall and Lions [11]). 

A widely applied class of numerical schemes for (1.1) is the class of finite difference schemes. In [12], 
Crandall and Lions proved the convergence of monotone finite difference schemes to the viscosity solutions of 

(1.1) . Unfortunately, monotone schemes are at most first order accurate, measured by local truncation errors 
in smooth regions of the solution. In Osher and Sethian [22] and Osher and Shu [23], a class of high order 
essentially non- oscillatory (ENO) schemes were introduced, which were adapted from the methodologies for 
hyperbolic conservation laws ([14], [25], [26]). The numerical results in [22] and [23] indicate convergence to 
the viscosity solutions of (1.1) with high order accuracy in smooth regions, and with sharp resolution of the 
discontinuous derivatives. Recently, weighted ENO (WENO) schemes for conservation laws ([20], [17]) have 
also been adapted to the Hamilton- Jacobi equations (1.1), ([15]). 

Finite difference methods require a structured mesh, hence are difficult to apply to complicated geometry 
or for adaptive mesh refinements. Finite volume schemes, which are based on arbitrary triangulation, are 
thus attractive in such cases. First order monotone type finite volume schemes were studied in [2]. A second 
order ENO type finite volume scheme is developed in [19]. However, higher order finite volume schemes face 
the problem of reconstruction on arbitrary triangulation, which is quite complicated. (See, e.g., [1]). 

The Runge-Kutta discontinuous Galerkin method ([5, 6, 7, 8, 9]) is a method devised to numerically 
solve the conservation law (CL): 

(1.2) Ut + fi(u) Xl -\ f fd(y)x d = o, u(x, 0) = u°(x). 
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The method has the following attractive properties: 

• It can be designed for any order of accuracy in space and time. In fact, p - version or spectral element 
type version can be designed ([21]); 

• It is an explicit method, thus efficient for solving the hyperbolic problem (1.2). No global lin ear or 
nonlinear systems must be solved; 

• It combines the flexibility of finite element methods in the easy handling of complicated geometry, 
with the high resolution property for discontinuous solutions of finite difference and finite volume 
methods through monotone fluxes or approximate Riemann solvers applied at the element interfaces 
and limiters; 

• It has nice stability properties: a local cell entropy inequality for the square entropy can be proven 
([16]) for general triangulation for any scalar nonlinear conservation laws (1.2) in any spatial dimen- 
sions and for any order of accuracy, without the need of nonlinear limiters. This implies nonlinear 
L 2 stability and entropy consistency even for discontinuous solutions; 

• The method is highly compact: the evolution of the information in any element depends only on 
the information of itself and its immediate neighbors, regardless of the order of accuracy. This is 
in contrast with high order finite volume schemes which must use wide stencils for the high order 
reconstruction. This compactness is responsible for the efficient parallel implementation of the 
method, see, e.g. [3], 

For details of the Runge-Kutta discontinuous Galerkin method, see the references listed above and also 
the review paper of Cockburn [4]. This method can also be generalized to solve problems containing higher 
derivatives ([10]), which is important for some of our applications in Sect. 4. 

It is well known that the Hamilton- Jacobi equation (1.1) is closely related to the conservation law (1.2), 
in fact in one space dimension d = 1 they are equivalent if one takes <j) = u x . It is thus not surprising that 
many successful numerical methods for the Hamilton- Jacobi equation (1.1) are adapted from those for the 
conservation law (1.2). The high order finite difference ENO methods in [22] and [23] are such examples. 

In this paper we adapt the Runge-Kutta discontinuous Galerkin method to solve the Hamilton- Jacobi 
equation (1.1). In Sect. 2 the algorithm in one space dimension is described and discussed. In Sect. 3 the 
algorithm for two space dimensions is developed. Numerical examples of one and two spatial dimensions are 
presented in Sect. 4. Concluding remarks are given in Sect. 5. 

2. One Dimensional Case. In one space dimension (1.1) becomes 

(2.1) </> t + H(<t> x ) = 0, *(*,O) = 0°(s). 

This is a relatively easy case because (2.1) is equivalent to the conservation law 

(2.2) u t -f H{u) x = 0, u(x,0) = u°(x) 


if we identify u = <f> x . 

Assuming (2.1) is solved in the interval a < x < b and it is divided into the following cells: 


(2.3) 

we denote 

(2.4) 


a = xi <i3 < 


< x 


AT+I 




= (*i - x i = \ ( x j-§ + x ]+h) > h i = x j+\ ~ h = maxhj, 
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and define the following approximation space 


(2.5) 




Here P k (Ij) is the set of all polynomials of degree at most k on the cell Ij. 

A k-th order discontinuous Galerkin scheme for the one dimensional Hamilton- Jacobi equation (2.1) can 
then be defined as follows: find ip G V k , such that 


(2.6) [ (p x tvdx - [ H(ip x ) 

Jij Jij 


v x dx + Ha , i v 




= 0, j = 


for all v e Here 


(2.7) 




is a monotone flux, i.e. H is non- decreasing in the first argument and non- increasing in the second, symbol- 
ically H( T,i), is Lipschitz continuous in both arguments, and is consistent, i.e. H(u,u ) = H(u). We will 
mainly use the simple (local) Lax- Friedrichs flux 

(2.8) H(u~,u + ) = H ^-^a(u+-M~) 

where a = max u ji/^u)! with the maximum taken over the range covered by u~ and For other monotone 
fluxes, e.g. the Godunov flux (see [23]). Notice that the method described above is exactly the discontinuous 
Galerkin method for the conservation law equation (2.2) satisfied by the derivative u — (j) x . This only 
determines (p for each element up to a constant, since it is only a scheme for <p x . The missing constant can 
be obtained in one of the following two ways: 

1. By requiring that 

(2.9) / (<p t + H (tp x )) vdx = 0, j = 1 , ..., N, 
for all v G V®, that is, 

(2.10) [ {ipt + H(<p x )) dx = 0, j = l,...,AL 

Jh 

2. By using (2.10) to update only one (or a few) elements, e.g., the left-most element 7i, then use 

(2.11) <p(xj,t) = p(xi,t) + / <p x (x,t)dx 

Jxi 

to determine the missing constant for the cell Ij. 

Both approaches are used in our numerical experiments. They perform similarly for smooth problems, with 
the first approach giving slightly better results. However, it is our numerical experience that, when there are 
singularities in the derivatives, the first approach will often produce dents and bumps when the integral path 
in time passes through the singularities at some earlier time. The philosophy of using the second approach 
is that one could update only a few elements whose time integral paths do not cross derivative singularities. 
The numerical results shown in Sect. 4 are obtained with the second approach. 

About the stability of the method proposed above, we can quote the following result of Jiang and Shu 
[16]. Here we assume compact support or periodic boundary condition for p. 
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Lemma 2.1. [16]. The following L 2 stability result for the derivative <p x holds for the discontinuous Galerkin 
method defined above, of any order of accuracy k applied to any nonlinear Hamilton- Jacobi equation (2.1): 

( 2 . 12 ) j t j\ 2 x dx< 0 . 

For a finite interval [a, 6], (2.12) trivially implies TVB (total variation bounded) property for the numer- 
ical solution <p\ 

(2.13) TV (ip) = jf \ip x \ dx < (^°( x )) dx - 

This is a rather strong stabihty result, considering that it applies even if the derivative of the solution 4> x 
develops discontinuities, no limiter has been added to the numerical scheme, and the scheme can be of 
arbitrary high order in accuracy. It also implies convergence of at least a subsequence of the numerical 
solution ip when h — ► 0. However, this stability result is not strong enough to imply that the limit solution 
is the viscosity solution of (2.1). 

Up to now we have only described the spatial discretization and have left the time variable t continuous 
(the method of fines approach). Time discretization is by the TVD (total variation diminishing) high order 
Runge-Kutta methods developed in [25] (see also [13]). The second and third order versions used in this 
paper are as follows: for solving the method of fines ODE 

(2.14) <p t = L(ip ), 
the second order TVD Runge-Kutta method is given by 

(2.15) v> (1) = p n + A tL(ip n ) 

<P n+1 = \<P n + \<P W + \a tL(<pW), 
and the third order TVD Runge-Kutta method is given by 

ip (l) = <p n + A tL(ip n ) 

(2.16) ^ 2 > = + ±AtL(<pW) 

¥> n+1 = + | v ( 2) + ^A tL(<pW). 

These Rimge-Kutta methods will also be used for the two dimensional case discussed in the next section. 

3. Two Dimensional Case. We consider in this section the case of two spatial dimensions. The 
algorithm in more spatial dimensions is similar. This time, the scalar Hamilton- Jacobi equation 

(3.1) <pt + H(4> x ,4> y ) = 0, tf>(x,y,0) = <j>°(x,y) 

is in some sense equivalent to the following conservation law system 

(3.2) u t + H(u,v) x = 0, v t + H(u,v) y = 0, (u, u)(x, y, 0) = (u, v)°(x, y). 

if we identify 

(3.3) («,v) = (<t> x ,<t>y). 
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For example, a vanishing viscosity solution of (3.1) corresponds, via (3.3), to a vanishing viscosity solution 
of (3.2), and vice versa ([18]). However, (3.2) is not a strictly hyperbolic system, which may cause problems 
in its numerical solution if we treat u and v as independent variables. Instead, we would like to still use (ft 
as our solution variable (a polynomial) and take its derivatives as u and v . 

Assuming we are solving (3.1) in the domain SI, which has a triangulation 7*. consisting of triangles or 
general polygons of maximum size (diameter) h, with the following approximation space 


(3.4) 


v£ = {v : v\k G P k (K), VATeTfc}, 


where P k (K) is again the set of all polynomials of degree at most k on the cell K. We propose a discontinuous 
Galerkin method for (3.1) as follows: find 93 € Vjj' , such that 


(3.5) 


and 

(3.6) 


/ <p xt vdxdy- / H(<p x , <p y )v x dxdy + V] / Hi, e ,K vdT = 0 

Jk Jk e edK Je 


/ (fytvdxdy- / H((p x ,ipy)vydxdy+ V' / H 2 ,e,Kvdr = 0 

Jk Jk e€dK J e 


for all v £ 1 and all K £ 7/ t , in a least square sense. Here the numerical flux is 

(3.7) = H i<etK ((V^) in<( *>, (Vv/r tW ) , t = 1,2 

where the superscript int(K) implies that the value is taken from within the element AT, and the superscript 
ext (K) implies that the value is taken from outside the element K and within the neighboring element K f 
sharing the edge e with K. The flux (3.7) satisfies the following properties: 

1. Hi, e ,K is Lipschitz continuous with respect to all its arguments; 

2. Consistency: 


JW(W,V0) = fr(V0K, 


where n = (ni, 712) is the unit outward normal to the edge e of the element K; 

3. Conservation: 

Hi,e,K ((V^r^MVyO^W) = -Hi,e,K’ (V93) e * t(ff,) ) . 

where K PI K f = e. 

We will again mainly use the simple (local) Lax- Friedrichs flux 

(3.8) H lte<K ((u,v)~ ,(u,v) + ) = H ni-^a(u + -u“) 

and 

(3-9) H 2 ,e,K((u, v)~, (u, v) + ) = H , - t -- ^ n 2 - ^/?(u + - v~) 

where a — max Uit; I 1 and (3 — max U)1J |, with the maximum being taken over the relevant 

(local) range. 

Notice that (3.5)- (3.6) are exactly the discontinuous Galerkin method for the conservation law system 
(3.2) satisfied by the derivatives (u,v) = V0. For a rectangular mesh and for k = 1 this recovers the (local) 
Lax-Friedrichs monotone scheme ([23]) for (3.1) if we identify Uij = anc ^ Vi i ~ ~ ’^ + 2Ay* ,J l 4 
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Of course, (3. 5)- (3. 6) have more equations than the number of degrees of freedoms (an over- determined 
system) for k > 1, thus a least square solution is needed. In practice, the least square procedure is performed 
as follows: we first evolve (3.2) for one time step (one inner stage for high order Runge-Kutta methods), 
using the discontinuous Galerkin method (3.5)- (3.6) with u = <p x and v = <p y \ then <p at the next time level 
(stage) is obtained (up to a constant) by least square: 


IKV’x - u) 2 + (<fy ~ v) 2 \\ L i(K) 


min 

^ eP k (K ) 


IIW’z - w) 2 + Wy 


- v) 2 \\l'(K)- 


This determines <p for each element up to a constant, since it is only a scheme for V<p. The missing constant 
can again be obtained in one of the following two ways: 

1. By requiring that 

(3.10) [ (<p t + H((p x ,tp y ))vdxdy = 0, 

Jk 

for all u G V® and for all K € 7/,, that is, 


( 3 - 11 ) 


L' 


(<p t + <p y )) dxdy = 0, 


VK eT h \ 


2. By using (3.11) to update only one (or a few) elements, e.g., the corner element(s), then use 


(3.12) 


t) ~ ip(A, t)+ f (<fix dx + (fi y dy ) 


to determine the missing constant. The path should be taken to avoid crossing a derivative discon- 
tinuity, if possible. 

In Sect. 4 we will only show numerical results obtained with the second approach, see the remarks in the 
previous section for the one dimensional case. 

We remark that the procedure discussed above is easily implemented in any triangulation, e.g., for both 
rectangles and triangles. 


4. Numerical Examples. Example 4.1. One dimensional Burgers 1 equation: 

(4.1) U, + <^= 0, — 1 < x < 1 

0) = — COS(7Tx) 

with periodic boundary conditions. 

The local Lax- Friedrichs flux (2.8) is used. At t — 0.5/7T 2 , the solution is still smooth. We list the errors 
and the numerical orders of accuracy in Table 4.1. We observe that, except for the P l case which seems to 
be only first order, P k for k > 1 seems to provide close to (k + l)-th order accuracy. 

At t = 3.5/V 2 , the solution has developed a discontinuous derivative. In Fig. 4.1, we show the sharp 
corner- like numerical solution with 41 elements obtained with P k for k = 1, 2, 3, 4. Here and below, the solid 
line is the exact solution, the circles are numerical solutions (only one point per element is drawn). 


Example 4.2. One dimensional equation with a non- convex flux: 


(4.2) 


< fit - COS (4> x + 1) = 0, -1 < X < 1 

0(a;, 0) = — cos(7rx) 


with periodic boundary conditions. 
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Table 4.1 

Accuracy for ID Burgers equation , t = 0.5/7 r 2 . 



p 1 

p 2 

p 3 

P i 

N 

L 1 error 

order 

L 1 error 

order 

L 1 error 

order 

L l error 

order 

10 

0.17E+00 

— 

0.14E-02 

— 

0.21E-03 

— 

0.57E-05 

— 

20 

0.78E-01 

1.12 

0.18E-03 

2.92 

0.13E-04 

3.94 

0.73E-06 

1 ^^ 

40 

0.35E-01 

1.16 

0.24E-04 

2.97 

0.75E-06 

4.17 

0.32E-07 


80 

0.16E-01 

1.12 

0.28E-05 

3.08 

0.43E-07 

4.12 

0.12E-08 


160 

0.76E-02 


0.31E-06 

3.19 

0.25E-08 

4.10 

0.48E-10 



p 1 

p 2 

P 3 

pi 

N 

L°° error 

order 

L°° error 


L°° error 

order 

L°° error 

order 

10 

0.29E+00 

— 

0.24E-02 

— 

0.69E-03 


0.13E-04 

— 

20 

0.13E+00 

1.13 

0.33E-03 


0.61E-04 

3.51 

0.16E-05 

2.99 

40 

0.58E-01 

1.15 

0.37E-04 

ESI 

0.58E-05 

3.39 

0.13E-06 


80 

0.27E-01 

1.11 

0.48E^05 

2.97 

0.38E-06 

3.93 

0.59E-08 

4.44 

160 

0.13E-01 

1.07 

0.59E-06 

3.00 

0.23E-07 

4.07 

0.25E-09 

4.57 


Table 4.2 

Accuracy for ID non-convex, H{u) = — cos(u + 1), t = 0.5/7T 2 . 



p 1 

p 2 

p 3 

P i 

N 

L l error 

order 

L l error 

order 

L 1 error 

order 

L 1 error 

order 

10 

0.84E-01 

- 

0.10E-02 


0.34E-03 

— 

0.24E-04 

— 

20 

0.36E-01 

1.23 

0.15E-03 

2.75 

0.30E-04 

3.49 

0.13E-05 

4.28 

40 

0.15E-01 



2.84 

0.15E-05 

4.33 

0.59E-07 

4.42 

80 

0.68E-02 




0.94E-07 

4.00 

0.21E-08 

4.78 


p 1 

p2 

P 3 

p4 

N 

L°° error 

order 

L°° error 

order 

L°° error 

order 

L°° error 

order 

10 

0.18E+00 

— 

0.15E-02 

— 

0.11E-02 

— 

0.99E-04 

— 

20 

0.73E-01 

1.31 

0.27E-03 


0.22E-03 

2.35 

0.13E-04 

2.95 

40 

0.31E-01 

1.24 


2.54 


3.63 

0.59E-06 


80 

0.14E-01 

1.16 

0.85E-05 

2.47 

0.14E-05 

3.75 

0.26E-07 

4.49 


The local Lax-Friedrichs flux (2.8) is used. At t — 0.5/7 r 2 , the solution is still smooth. The accuracy of 
the numerical solution is listed in Table 4.2. We observe similar accuracy as in the previous example. 

At t = 1.5/ 7 r 2 } the solution has developed corner- like discontinuity in the derivative. The numerical 
result with 41 elements is shown in Fig. 4.2. 

Example 4,3. Riemann problem for the one dimensional equation with a non-convex flux: 

, 43) { + l)(£-4)=0, — 1 < x < 1 

\ 0(*,O) = -2|*| 

For this test problem, the discontinuous Galerkin method, as described in Sect. 2 and applied in the 
previous two examples, fails to “open up” the initial single discontinuity in the derivative sufficiently to 
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Fig. 4.1. One- dimensional Burgers’ equation, t — 3.5/n 2 . 


generate the correct entropy solution. We have found out that a nonlinear total variation bounded limiting, 
described in detail in [6] for this one dimensional case, is needed for convergence towards the entropy solution. 
This and the 2D Riemann problem in Example 4.6 below are the only two examples in this paper in which 
we use the nonlinear limiting. We remark that for the finite difference schemes, such nonlinear limiting or 
the adaptive stencil in ENO is needed in most cases in order to enforce stability and to obtain non-oscillatory 
results. 

Numerical results at t — 1 with 81 elements, using the local Lax-Friedrichs flux (2.8), is shown in Fig. 
4.3. The results of using the Godunov flux is shown in Fig. 4.4. We can see that while for P 1 , the results of 
using two different monotone fluxes are significantly different in resolution, this difference is greatly reduced 
for higher order of accuracy. In most of the high order cases, the simple local Lax-Friedrichs flux (2.8) should 
be good enough. 


Example 4.4. Two dimensional Burgers’ equation: 


(4.4) 


4>t + ( - 0 * + * 2 v+l)2 =0, -2 < x < 2, -2 < y < 2 

4>(x, y, 0) = - cos j 


with periodic boundary conditions. 
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Fig. 4.2. One dimension non-convex, H(u) = — cos(ti + 1), t = 1.5/7T 2 . 

We first use uniform rectangular meshes and the local Lax- Friedrichs flux (3.8)- (3.9). At t — 0.5/7T 2 , the 
solution is still smooth. The errors and orders of accuracy are listed in Table 4.3. It seems that only k- th 
order of accuracy is achieved when (p is a piecewise polynomial of degree k. 

At t = 1.5/n 2 , the solution has discontinuous derivatives. Fig. 4.5 is the graph of the numerical solution 
with 40 x 40 elements. 

Next we use triangle based triangulation, the mesh with h = \ is shown in Fig. 4.6. The accuracy at 
t = 0.5/7T 2 is shown in Table 4.4. Similar accuracy pattern is observed as in the rectangular case. The result 
at t = 1.5/7T 2 , when the derivative is discontinuous, is shown in Fig. 4.7. 

Example 4.5. Two dimensional equation: 

{ 4>t — cos(<Ae -j- <j>y + 1) = 0, “2 < x < 2, —2 < y < 2 

^ 4 ' 5 ^ \ 4>{x,y,Q) = -cos(^b^) 

with periodic boundary conditions. 

For this example we use uniform rectangular meshes. The local Lax- Friedrichs flux (3.8)- (3.9) is used. 
The solution is smooth at t — 0.5/7 r 2 . The accuracy of the numerical solution is shown in Table 4.5. 

The solution has developed a discontinuous derivative at t = 1.5/7T 2 . Results with 40 x 40 elements are 
shown in Fig. 4.8. 
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Table 4.3 

Accuracy for 2D Burgers equation , rectangular mesh , t = 0.5/7T 2 . 



p 1 

P 2 

P 3 

TV x TV 

L 1 error 

order 

L 1 error 

order 

L 1 error 

order 

10 x 10 



8.62E-03 

— 

3.19E-03 

— 

20 x 20 

3.36E-02 

1.268 

1.72E-03 

2.325 

3.49E-04 

3.192 

40 x 40 

1.48E-02 

1.183 

3.93E-04 

2.130 

6.64E-05 

2.394 

80 x 80 

6.88E-03 

1.105 

9.74E-05 

2.013 

1.14E-05 

2.542 

160 x 160 

3.31E-03 

1.056 

2.45E-05 

1.991 

1.68E-06 

2.763 


P 1 

p 2 

P 3 

TV x TV 

L°° error 

order 

L°° error 

order 

L°° error 

order 

10 x 10 

2.62E-01 


3.56E-02 

— 

8.65E-03 

— 

20 x 20 

1.14E-01 


8.40E-03 

2.083 

1.16E-03 

2.899 

40 x 40 

5.00E-02 


2.02E-03 

2.056 

1.98E-04 

2.551 




4.92E-04 

2.038 

3.13E-05 

2.661 

160 x 160 

1.16E-02 

1.043 

1.21E-04 

2.024 

4.41E-06 

2.827 


Table 4.4 

Accuracy for 2D Burgers equation, triangular mesh, t = 0.5/n 2 . 



P 2 

P 3 

h 

L l error 

order 

L°° error 

order 

L l error 

order 

L°° error 

order 

1 

5.48E-02 

— 

1.52E-01 

— 

1.17E-02 

--- 

2.25E-02 

— 

1/2 

1.35E-02 


6.26E-02 

1.28 

1.35E-03 

3.12 

4.12E-03 


1/4 

2.94E-03 


1.55E-02 

2.01 

1.45E-04 

3.22 

4.31E-04 


1/8 

6.68E-04 

2.14 

3.44E-03 

2.17 

1.71E-05 

3.08 

7.53E-05 

2.52 


Table 4.5 

Accuracy, 2D, H(u , u) = — cos(u + v + 1), t = 0.5/7T 2 . 



p 1 

P 2 

P 3 

TV x TV 

L 1 error 

order 

L 1 error 

order 

L 1 error 

order 

10 x 10 

6.47E-02 

— 

8.31E-03 

— 

1.35E-02 

— 

20 x 20 

2.54E-02 

1.349 

1.93E-03 

2.106 

1.57E-03 

3.104 

40 x 40 

1.05E-02 

1.274 

4.58E-04 

2.075 

2.39E-04 

2.716 

80 x 80 

4.74E-03 

1.147 

1.13E-04 

2.019 

2.89E>05 

3.048 

160 x 160 

2.23E-03 

1.088 

2.83E-05 

1.997 

4.38E-06 

2.722 


p 1 

p 2 

P 3 

TV x N 

L°° error 

order 

L°° error 

order 

L°° error 

order 

10 x 10 

1.47E-01 

— 

1.88E-02 


2.36E-02 


20 x 20 

6.75E-02 

1.123 

7.34E-03 

1.357 

3.44E-03 

2.778 

40 x 40 

2.65E-02 

1.349 

1.83E-03 

2.004 

4.59E-04 

2.906 

80 x 80 

1.18E-02 

1.167 

4.55E-04 

2.008 

5.78E-05 

2.989 

160 x 160 

2.23E-03 

1.088 

1.13E-04 

2.010 

8.54E-06 

2.759 
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Fig. 4.3. One dimension Riemann problem, local Lax- Friedrichs flux , H(u) — |(u 2 — l)(u 2 — 4), f — 1. 


Example 4.6. Two dimensional Riemann problem: 


(4.6) 


<f>t + sin(0 x + <t> y ) = 0, -1 < x < 1, -1 < y < 1 

(f>(x,y, 0 ) = tt (|?/| - |*|) 


For this example we use a uniform rectangular mesh with 40 x 40 elements. The local Lax- Friedrichs 
flux (3.8)- (3.9) is used. As it was mentioned in Example 4.3, we have found out that a nonlinear limiting 
[8], [9] is needed, for convergence towards an entropy solution. We show the numerical solution at t = 1 in 
Fig. 4.9. 


Example 4.7. The problem of a propagating surface: 


(4.7) 


4>t - (1 - eK) yjl + $1 + ^ =0, 0<x<l,0<i/<l 

4>{x, y, 0) = 1 - |(cos(27rx - 1)) (cos( 27 tj/ - 1)) 


where K is the mean curvature defined by 

< ^ > xx(l ~b (fry) ~~ 2(f) X y(j) x (f)y + (f)yy{l ~b <f> x ) 


(4.8) 


K = 


(1 + £ + **)* 

and e is a small constant. Periodic boundary condition is used. 


ll 




Fig. 4.4. One dimension Riemann problem, Godunov flux, H(u) — \(u 2 — l)(u 2 —4 ),t = 1. 


This problem was studied in [22] by using the finite difference ENO schemes. We apply the discontinuous 
Galerkin method, with the second derivative terms handled by the local discontinuous Galerkin techniques 
presented and analyzed in [10], which amounts to solving the following system 

' u t - (vT + u 2 + ^ 2 + e p(1+ " VuC* r(1+u2) ) x = 0 
Vt - (vi + ^ + t ; 2 + £ p( 1 + " 2 \+ 2 u q O r(1+tt3) ) 1 = o 

( 4 ‘ 9 ) p- Ul = 0 

q - Uy = 0 

r — v y = 0 

using the discontinuous Galerkin method. The details of the method, especially the choices of fluxes, which 
are important for stability, can be found in [10]. 

We first use a uniform rectangular mesh of 50 x 50 elements and the local Lax- Friedrichs flux (3. 8)- (3.9). 
The results of e = 0 (pure convection) and € = 0.1 are presented in Fig. 4.10 and Fig. 4.11, respectively. 
Notice that the surface at t = 0 is shifted downward by 0.35 in order to show the detail of the solution at 
t = 0.3. 

Next we use a triangulation shown in Fig. 4.12. We refine the mesh around the center of domain where 
the solution develops discontinuous derivatives (for the € = 0 case). There are 2146 triangles and 1108 nodes 
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P , 40x40 elements 



P , 40x40 elements 



Fig. 4.5. Two dimension Burgers’ equation , rectangular mesh , t=1.5/n 2 . 



FlG. 4.6. Triangulation for two dimensional Burgers equation, h = 


in this triangulation. The solutions are displayed in Fig. 4.13 and Fig. 4.14, respectively, for e = 0 (pure 
convection) and e = 0.1. Notice that we again shift the solution at t = 0.0 downward by 0.35 to show the 
detail of the solutions at later time. 


Example 4.8. The problem of a propagating surface on a unit disk. The equation is the same as (4.7) in 
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P 2 , h = 1/8 


P 3 , h = 1/8 



Fig. 4.7. Two dimension Burgers’ equation, triangular mesh., t=1.5/ n 2 . 



Fig. 4.8. Two dimensional , H(u, v) = — cos(u + v + 1), t = 1.5/7T 2 . 

the previous example, but it is solved on a unit disk x 2 + y 2 < 1 with an initial condition 

/ 7 t(x 2 + y 2 ) \ 

0(x,y,O) = sin ( J 

and a Neumann type boundary condition V<^ = 0. 

It is difficult to use rectangular meshes for this problem. Instead we use the triangulation shown in Fig. 
4.15. Notice that we have again refined the mesh near the center of the domain where the solution develops 
discontinuous derivatives. There are 1792 triangles and 922 nodes in this triangulation. The solutions with 
£ = 0 are displayed in Fig. 4.16. Notice that the solution at t = 0 is shifted downward by 0.2 to show the 
detail of the solution at later time. 
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Fig. 4.10. Propag 


solution with e = 0.1 are displayed 
d by 0.2 to show the detail of the s< 





P a , 50x50 elements 


P 3 , 50x50 elements 



Fig. 4.11. Propagating surfaces, rectangular mesh , e = 0.1. 


Example 4.9. A problem from optimal control [23]: 



4>t + (sin y)<t> x + (sinx + sign (<j> y ))<i>y - \ sin 2 y - (1 - cosx) = 0, -n < x < 7r, -n < y < tt 

4>{x,V, 0) = 0 


with periodic boundary conditions. We use a uniform rectangular mesh of 40 x 40 elements and the local 
Lax-Friedrichs flux (3.8)-(3.9). The solution at t = 1 is shown in Fig. 4.18, while the optimal control 
w = sig n(4> y ) is shown in Fig. 4.19. 

Notice that our method computes V0 as an independent variable. It is very desirable for those problems 
in which the most interesting features are contained in the first derivatives of 0, as in this optimal control 
problem. 


Example 4.10. A problem from computer vision [24]: 

f <Pt + I(x,y)yJ\ + ^ + (^2 -1 = 0, -1 < X < 1, -1 < y < 1 

1 i/,0) = 0 

with 0 — 0 as the boundary condition. The steady state solution of this problem is the shape lighted by 
a source located at infinity with vertical direction. The solution is not unique if there are points at which 
I(x, y) — 1. Conditions must be prescribed at those points where I(x,y) = 1. Since our method is a finite 
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Fig. 4.12. Triangulation used for the propagating surfaces. 

element method, we need to prescribe suitable conditions at the correspondent elements. We take 

(4.12) I(x, y) = 1/Vl + (1-M) 2 + (1-M) 2 

The exact steady solution is </>(x,y, oo) = (1 — |x|)(l — |y|). We use a uniform rectangular mesh of 40 x 40 
elements and the local Lax- Friedrichs flux (3.8)-(3.9). We impose the exact boundary conditions for u = 
0 X , v = 4> y from the above exact steady solution, and take the exact value at one point (the lower left corner) 
to recover cf>. The results for P 2 and P 3 are presented in Fig. 4.20, while Fig. 4.21 contains the history of 
iterations to the steady state. 

Next we take 

(4.13) 7(x,y) = l/y/l + 4y 2 (l - x 2 ) 2 + 4x 2 (l - y 2 ) 2 

The exact steady solution is 0(x, y, oo) = (1 — x 2 ) (1 — y 2 ). We again use a uniform rectangular mesh of 40 x 40 
elements, the local Lax- Friedrichs flux (3.8)-(3.9), impose the exact boundary conditions for u = <f> xi v = <j) y 
from the above exact steady solution, and take the exact value at one point (the lower left corner) to recover 
<j>. A continuation method is used, with the steady solution using 

(4.14) 7 e (x,y) = l/y/l + 4y 2 (l - x 2 ) 2 -F 4x 2 (l - y 2 ) 2 + e 

for bigger e as the initial condition for smaller e . The sequence of e used are e = 0.2,0.05, 0. The results for 
P 2 and P 3 are presented in Fig. 4.22. 

5. Concluding Remarks.. We have developed and tested a class of discontinuous Galerkin methods 
for solving nonlinear Hamilton- Jacobi equations. These methods have the advantage of easy handling of 
complicated geometry and local mesh refinements, as well as efficient parallel implementations. Numerical 
examples in one and two space dimensions are shown to illustrate the capability of the methods. 
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Fig. 4.13. Propagating surfaces, triangular mesh, e = 0. 
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P 2 , 40x40 elements 


P 3 , 40x40 elements 




Fig. 4.20. Computer vision problem, <f>(x,y, oo) = (1 — |x|)(l — |y|). 




Fig. 4.21. Computer vision problem, history of iterations. 
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40x40 elements 


P 3 , 40x40 elements 



Fig. 4.22. Computer vision problem , 4>(x, y, oo) = (1 -i 2 )(l - y 2 ). 
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